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Abstract 

Context. Active Galactic Nuclei are highly variable emitters of electromagnetic waves from the radio to the gamma-ray regime. This 
variability may be periodic, which in turn could be the signature of a binary black hole. Systems of black holes are strong emitters 
of gravitational waves whose amplitude depends on the binary orbital parameters as the component mass, the orbital semi-major-axis 
and eccentricity. 

Aims. It is our aim to prove the existence of periodicity of the AGN Markarian 501 from several observations in different wavelengths. 
A simultaneous periodicity in different wavelengths provides evidence for bound binary black holes in the core of AGN. 
Methods. Existing data sets from observations by Whipple, SWIFT, RXTE and MAGIC have been analysed with the Lomb-Scargle 
method, the epoch folding technique and the SigSpec software. 

Results. Our analysis shows a 72-day period, which could not be seen in previous works due to the limited length of observations. 
This does not contradict a 23-day period which can be derived as a higher harmonic from the 72-day period. 
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1. Introduction 

Active Galactic Nuclei (AGN) have been a class of most inter- 
esting astronomical sources for a long time, mainly for their 
high variability on time scales from years down to minutes. 
Nowadays these variations may be observed from radio up to 
gamma-rays (> 10 TeV). In most cases the variation of the AGN 
do not show signs of periodicity but a unpredictable switch- 
ing between quiescent and active states and flaring behaviour 
The explanation of such behaviour requires several different ap- 
proaches in order to explain the very different timescales. Most 
of the models make use of the source size (the black hole radius) 
for the minimum time scales or the accretion efficiency and outer 
scales for the maximum time scales. 

Some sources show a periodicity on a timescale of several 
tens of days in the optical. X-ray or TeV data like Mrk 421, 
Mrk 501, 3C66A and PKS 2155-304 JHayashida et al.||1998[ 



ILainela et al.[T999l|Kranich et al.|2001|r6sone & Teshima|200I| l, 
whereas the long-term optical light curves from classical sources 
such as BL Lacertae, ON 231, 3C273, OJ 287, PKS 0735-H178, 
3C345 and AO 0235-1-16 usually suggest (observed) timescales 
of several years (e.g. |Sillanpaa et al.|[T988} [Fan et aL]|I997 



|Raiteri et al.|200l} . Models invoked to describe aperiodic vari- 
ability do not give rise to explanations for periodic variations. It 
is therefore necessary to develop new models. One of these mod- 
els is the binary black hole (BBH) model, which is explained in 
detail in the next section. It is worth noting that the origin of 
these quasi periodic objects (QPOs) may be related to instabili- 
ties of (optically thin) accretion disks (cf. e.g. [Fan e t al. (2008jl 
and citations within). We do not want to stress this point but con- 
centrate on the BBH model to explain the variability. 

The existence of BBHs seems plausible since the host galax- 
ies of most AGN are elliptical galaxies which originate from 
galaxy mergers. The study of BBHs and their interactions be- 



come of primary importance from both the astrophysical and 
theoretical points of view now that the first gravitational wave 
detectors have started operating; gravitational waves, while cer- 
tainly extremely weak, are also extremely pervasive and thus a 
novel and very promising way of exploring the universe. For a 
review on BBHs the reader is referred to Komossa (2003). 
BBHs have already been detected ([Hudson et al. 



2006 



Sudou|i2003| |Sillanpaa|[T998t [Romero et al.||2003[ ), but most 
sources are far too distant to be separated in radio or x-rays 
(for a detailed discussion of the emission mechanisms see 
(Bogdanovic et al. 2008| l. Using a multi-messenger technique 
(combined electromagnetic and gravitational wave astronomy) 
it should be possible to identify periodic AGN as BBH systems 
by means of their gravitational wave signal on the one hand 
while on the other hand x-ray and gamma-astronomy could eas- 
ily identify possible targets for gravitational wave telescopes. 

The aim of this paper is to present a complete data analysis of 
Markarian 501 in different wavelengths searching for periodic- 
ity. This will allow for restricting possible black hole masses and 
the distance between the two black holes. Taking into account 
the model of Rieger & Mannheim) f2000), where the smaller 
companion is emitting a jet with Lorentz factor y towards the 
observer, we expect periodic behaviour over the complete AGN 
spectrum. 

This paper is organised as follows. In Section|2]we review a 
model used to explain the variability and present the equations 
used to determine the parameters of the binary system under 
study. Then in Section [3] the data we used and the experiments 
are described. In Section |4] we describe methodology for data 
analysis and in|5]we present the results of the analysis, discuss 
our findings and draw some conclusions. 

2. The binary black hole model 
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As discussed above, the presence of close supermassive bi- 
nary black hole (SMBBH) systems has been repeatedly invoked 
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as plausible source for a number of observational findings in 
blazar-type AGN, ranging from misalignment and precession of 
jets to helical trajectories and quasi-periodic variability (cf. e.g. 
[Kom ossa 2006 ). 

What makes the BBH model an unique interpretation is that 
it facilitates facts othei-wise not possible: (1.) It is based on quite 
general arguments for bottom-up structure formation, (2.) it in- 
corporates helical jet trajectories observed in many sources, (3.) 
it provides a reasonable explanation for long term periodic vari- 
abiUty, (4.) it can explain, to some extend, quasi-periodic vari- 
abiUty on different time-scales in different energy bands (e.g. 



[Rieger 2007). 

The nearby TeV blazar Mr-k 501 (z = 0.033) attracted atten- 
tion in 1997, when the source underwent a phase of high activity 
becoming the brightest source in the sky at TeV energies. Recent 
contributions like Rieger & Mannheim (2000, 2001) have pro- 
posed a possible periodicity of Pobs - 23 days, during the 1997 
high state of Mrk 501, which could be related to the orbital mo- 
tion in a SMBBH system. The proposed model relies on the fol- 
lowing assumptions: (1.) The observed periodicity is associated 
with a relativistically moving feature in the jet (e.g. knot). (2.) 
Due to the orbital motion of the jet-emitting (secondary) BH, 
the trajectory of the outward moving knot is expected to be a 
long drawn helix. The orbital motion is thus regarded as the un- 
derlying driving power for periodicity. (3.) The jet, which dom- 
inates the observed emission, is formed by the less massive BH. 
(4.) The observed periodicity arises due to a Doppler-shifted flux 
modulation which is caused by a slight change of the inclination 
angle with respect to the line of sight. 

Following Rieger & Mannheim (2000)(hereinafter RM), we 
assume that the observed signal periodicity has a geometrical 
origin, being a consequence of a Doppler-shifted modulation. It 
is therefore possible to relate the observed signal period Poi,, to 
the Keplerian orbital period 



VG(m + M) 

^372 ' 



(1) 



by the equation (cf . |Camenzind & Krockenberger|1992[ |Roland| 
let al.|I994| ) 

Pohs = (l+Z){l-^COSiyK, (2) 

where v.. is the typical jet velocity ass umed to be - 



c(l - Jf, ) . Current emission models (e.g. 
an inclination angle ; =s 1 jjb , with typical bulk Lorentz 
in the range 10 
[1999). 

For a resolved emission region (e.g. a blob of plasma) with 
spectral index a, the spectral flux modulation by Doppler boost- 
ing can be written as 



Spada 



1999|) favour 
factors 



15 (e.g. |Mannheim e~al. 1996; Spada et al. 



(3) 



where S ' is the spectral flux density measured in the co-moving 
frame and 6{t) denotes the Doppler factor given by 



1 



rJi-^Si cos 0(f)]' 



(4) 



where 0{t) is the angle between the direction of the emission re- 
gion (e.g. blob) and the line of sight. Due to the orbital motion 
around the centre-of-mass, the Doppler factor for the emission 
region is a periodical function of time. In the simplest case (with- 
out determination) the Doppler factor may be written as follows 



(RM) 



^1-(v2 + Q2/;2)/c2 

1 - (Vj COS / + Q.iiR sin / sin ilj.?) /c 



(5) 



A periodically changing viewing angle may thus naturally lead 
to a periodicity in the observed lightcurves even for an intrinsi- 
cally constant flux. 

Depending on the position of the jet-emitting BH along its 
orbit, the Doppler factor has two extremal values, so that one 
obtains the condition 



l/(3+ff) 



(6) 



where / = S maySy) I S mm{v) is the observed maximum to min- 
imum ratio. Consequently, by inserting Eq. (|5]l in Eq. (j6]l one 
finds 
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Finally, by using the definition of the Keplerian orbital period 
Eq. ([T]) and noting that R - Mdj (m + M) one derives the binary 
mass ratio equation (RM) 



M 



(m + M)2/3 [(l+z)27rG]'/V'^^'^"' 
and from Eq. Q we get 
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To be able to determine the primary and secondary black hole 
masses m and M by solving simultaneously Eq. ([8]) and (j9]l 
we have to provide an additional assumption or a third equa- 
tion. RM obtained the binary separation d by equating the time- 
scale for gravitation radiation with the gas dynamical time-scale 
(Begelman et al. 1980 1. We have not adopted this condition 
since it is somewhat arbitrary to assume Tg^s - T(jw, nowadays. 
Therefore we will leave the distance d asa free parameter in the 
equations. Instead, we solve Eq. (|8]) and (j9]l by adopting reason- 
able values for the total binary mass frorn obser vational findings, 
M + m= 10^-^^ solar masses ( |Woo et al.|2005| ). 

As outlined, we can determine the orbital parameters of the 
system in case of a perfectly circular orbit. Notably, the assump- 
tion of a perfectly circular orbit is simplifying the problem to 
some extend. It is true that the orbits tend to circularise due to 
gravitational radiation, but this happens within a time-scale of 
the same order of magnitude as the merging time-scale (Peters] 
& Mathews|1963|[Fitchett 1987) . Therefore it is possible that the 
constituting black holes of a binary system are still on eccentric 
orbits, in which case the method by RM will be only an approx- 
imation to the real situation. Recently |De Paolis et al.| ( |2002| 
2003 ) extended the model of RM to elliptical orbits including 
the eccentricity e, which reduces to the circular model in case of 
e ^ 0. In case of an eccentric orbit a degeneracy in the solutions 
appears making it impossible to determine the orbital parameters 



in an unique fashion. As pointed out further by De Paolis et al 



(2002 ) the degeneracy can only be solved by the detection of a 
gravitational wave signal. 
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Excerpt from RXTE 50901-54803 ASM data 

Sample MJD 54000-54150 



SWIFT: 53414-54721 data from the BAT instrument 

Lightcurvi;MarkLirian50l 




Figure 1. Subset of the whole available ASM data on Mrk 501. 
The sampling is not equidistant and the data points clearly differ 
significantly from zero. 



3. Data 

Data was taken from a variety of different sources, earthbound 
telescopes Uke MAGIC or WHIPPLE, but also from satellite 
experiments like RXTE and SWIFT. In contrast to any previ- 
ous studies on this source, we explicitly excluded the time of 
the flare in 1997, so to determine whether a periodicity would 
still be present in the quiescent state of Mrk 501. We didn't 
make use of any raw data material, but only used data that was 
background-corrected and processed by the respective collabo- 
rations , to which we refer for more details. Therefore the units 
of the plots are arbitrary, a fact that doesn't cause any problems 
for our analysis as we are not interested in any absolute flux val- 
ues of the source. 



3.1. RXTE 

The Rossi X-ray Timing Explorer (RXTE) is an orbital X-ray 
telescope for the observation of black holes, neutron stars, pul- 
sars and X-ray binaries. The RXTE has three instruments on 
board: the Proportional Counter Array (PCA, 2-60 keV), the 
High-Energy X-ray Timing Experiment (HEXTE, 15 - 250 keV) 
and the so-called All Sky Monitor (ASM, 2 - 10 keV). The PCA 
is an array of five proportional counters with a time resolution of 
Ijus and an angular of 1° FWHM, the energy resolution amounts 
to 18% at 6 keV. For HEXTE the time sampling is 8jus with a 
field of view of 1° FWHM at an energy resolution of 15% at 60 
keV and the ASM consists of three wide-angle cameras which 
scan 80% of the sky every 90 minutes ^ahoda et al.|1996[[RXTE 
|2009i ). 

The data from RXTE we used was gathered by the ASM 
from March 1998 until December 2008 and a sample from the 
whole data set can be seen in Fig.[T] We are only plotting a subset 
to show the actual sampling structure of the RXTE data, firstly 
because the data is extremely plentiful (it contains 3697 data 
points) and secondly to have a direct comparison to the SWIFT 
data in Fig. [5] The data shows a rms of 5.15 ■ 10"' (in arbitrary 
units) with an average signal-to-noise ratio of 1.22. 




53600 53800 



54000 54200 54400 
MJD 



Figure!. Complete available data for SWIFT containing 1153 
data points with error bars. Background has already been sub- 
tracted before data retrieval (for details see (Swift 2009a) ). 



3.2. SWIFT 

The SWIFT satelhte is part of the SWIFT Gamma-Ray Burst 
Mission planned and executed by NASA dedicated to the study 
of GRBs and their afterglows. It consists of three instruments, 
the Burst Alert Telescope (BAT) which detects gamma-rays in a 
range from 15 to 150 keV with a field of view of about 2 sr and 
computes their position on the sky with arc -minute positional 
accuracy. After the detection of a GRB by the BAT, its image and 
spectrum can be observed with the X-Ray Telescope (XRT, 300 
eV - 10 keV) and the Ultraviolet/Optical Telescope (UVOT, 170 



- 650 nm) with an arc-second positional accuracy ( Barthelmy & 
ISwift Team|2003i|Swift|20"09bl l. 

The data that was used can be obtained from a database at the 
website of the NASA ( Swift 200 8). We made use of data from 
February 2005 to September 2008. As can be seen in Fig. [2] the 
lightcurve provided by the SWIFT satellite looks very much like 
background noise. Especially when comparing the size of the 
error-bars to the amplitude of the signal in the subset in Fig. [3] 
common sense implies that this lightcurve should be regarded 
with caution. The subset was chosen around the maximum am- 
plitude measured within the whole sample and the same time- 
span was then plotted for RXTE in Fig. [T]for comparison. The 
rms of the flux of the whole data sample is 35.7 ■ 10""* (in arbi- 
trary units) with an average signal-to-noise ratio of only 0.570. 
As a result of this, we consider the data from SWIFT merely as 
an indicator, but not as solid evidence that seriously influences 
the outcome of the hypothesis tests. 

3.3. MAGIC 

MAGIC is an Imaging Air Cherenkov Telescope (lACT) sit- 
uated at the Roque de los Muchachos on one of the Canary 
Islands, La Palma. It has an active mirror surface of 236 m^ and 
a hexagonal camera with a diameter of 1.05 m consisting of 576 
photomultipliers. The trigger and analysis threshold are 50 - 70 
GeV and the upper detection limit is 30 TeV, for further details 
on the MAGIC teles cope see |Petry & The MAGIC Telescope 
|Collaboration ( 1999| l. The data we used was collected from May 
to July 2005 in the energy range of 150 GeV to 10 TeV ( |Albert| 
et al. 20 07) and can be seen in Fig. |4] It contains of 24 data 
points with a rms of 2.63 ■ 10"'" (in arbitrary units) and an aver- 
age signal-to-noise ratio of 1 1 .0. 
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Figure 3. Subset of SWIFT data around the maximum ampli- Figure 5. Preliminary High-Energy data measured by 
tude in order to show the single points more clearly. The highest WHIPPLE from 2006 to 2008, the set contains 33 data 
peak lies 12.8 a above the mean flux and is therefore considered points 
a fluctuation. So the remaining of the lightcurve is believed to 
consist merely of zeroes. 
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Figure 4. Measured High-Energy data from MAGIC in 2005, the 
set contains 24 data points. 



3.4. WHIPPLE 

In the early 1980s, the 10 m Whipple Telescope was built at the 
Fred Lawrence Whipple Observatory in southern Arizona. It is 
still fully operational and is nowadays used for long-term mon- 
itoring of some interesting astrophysical sources in the energy 
range from 100 Gev to 10 TeV. The data of MrkSOl we used 
was taken from M arch 2006 to March 2008 and is preliminary 
( |WHIPPLE| |2008|. It shows a rms of 2.33 • 10~' (in ai-bitrai-y 
units) and an average signal-to-noise ratio of 2.45 (see Fig.|5]l. 
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4. Analysis technique 

In observational sciences the obsei-ver can often not completely 
control the time of the observations, but must simply accept a 
certain dictated set of observation times f,. In order to process un- 
evenly sampled data, simple Fourier-transformations cannot be 
used, but must be replaced by a more elaborate analysis. There 
are several ways of how to tackle this problem. In our work we 
made use of three different methods. These we selected by the 
criterion of being mutually fairly orthogonal so as to reduce sys- 
tematic errors. In the following paragraphs two of them will be 
outlined briefly, the third was gratefully accepted from the Dept. 
for Astroseismology of the University of Vienna. 

4. 1 . Lomb-Scargle and Epoch Folding 
4.1.1. Lomb-Scargle 

A full derivation of the Lomb normalised periodogram (spectral 
power as a function of angular frequency w = 2nf > 0) may 
be found in |Lornb) ([1976J and |Scargle| (| 1982) . The algorithm 
transforms a discretely and unevenly sampled set of data points 
into a quasi-continuous function in Fourier space, essentially by 
applying a ;t'^-test on an equivalent to linear least square fitting 
to the model 



h(t) - A cos (jjt + B sin ut 



(10) 



In order to detect low frequency signals, the oversampling 
method ( [Press et aL 2003| l is usually considered. However, spuri- 
ous effects are thus introduced and particular attention is neces- 
sary in analysing the resulting periodogram (see also Vaughan 
|2005 V An ill-chosen value of the oversampling parameter is 
found to produce artificial oscillations in the periodogram that 
can easily be interpreted wrongly. We used the oversampling 



method and in Sec. 4.3 describe means of controlling these os- 
cillations and we depict a way of statistically evaluating them as 
one of the major contributions to the uncertainty of an identified 
period. 

4.1 .2. Epoch Folding 



Following [Leahy et al.| ( |1983) , the epoch-folding statistic is de- 
fined as 



M 



(11) 



where M denotes the number of phase bins and cr^ the pop- 
ulation variance of Xj. A high value (» M - 1) will sig- 
nal the presence of a periodicity for which t he significance 
can be estimated from the x\K-^ distribution (cf. [Larsson|l996 



Schwarzenberg-Czerny 1 1 997| . Therefore, this method yields an 



estimate of whether or not the data is distributed Gaussianly; in 
the case of coloured noise being part of the experimental error, 
as is the case in the present study, the method may not yield re- 
liable results, as the noise itself, if the signal to noise ratio is too 
bad, could cause the distribution of the data to be non-Gaussian. 
We used this method especially in order to cross-check the re- 
sults of the other algorithms for artifacts. This means, whenever 
a periodogram shows high values in Fourier power, we construct 
the phase diagram of the corresponding period and try fits onto 
formulae of periodic signals, like: 



y - Aq+ Al ■ sin(A2 ■ x + A3) 



(12) 



where the A,- are the degrees of freedom in the fit. We also tried 
higher derivatives or higher powers of trigonometric functions, 
however the formula Eq.[T2[ proved to yield the best results. 

4.2. SigSpec 

This software (for details see |Reegen[[2007| l presents a combi- 
nation of the techniques mentioned above, and has many more 
other features in addition. One aspect that needs to be mentioned 
is that where the Lomb-Scargle Periodogram is presented as a 
continuous function of the period, SigSpec only provides dis- 
crete values of those periods, that are above a certain threshold- 
level of spectral significance. The latter is the inverse False- 
Alarm-Probability f(> 37) = l-(l-exp(-y))*' scaled logarithmi- 
cally (specified significance threshold y and number of indepen- 
dent frequencies M). The conversion of a chosen threshold for 
maximum spectral significance into its corresponding "individ- 
ual" spectral significance threshold for accepting detected peaks 
is very similar to |Scargle (1982) and can be found in Reegen] 
( 2007[ l. SigSpec cannot deal with coloured noise, however it is 
possible to choose different spectral significance thresholds for 
different frequency regions (Reegen_2007 ). 



4.3. Random-noise-modelling and Random-sampling 

In order to handle the uncertainties of the results we used two ap- 
proaches as follows. One way of testing the significance of the 
resulting Fourier Lomb-Scargle-transformation is "Random- 
sampling": we randomly choose a certain number P of data 
points from the original data set thus creating n subsets, trans- 
form all of them (using one of the three algorithms at a time) 
individually and add the thus acquired n Fourier-sets together 
so to form one single periodogram. P was chosen to be 2/3 of 
the original length of the lightcurve so as to guarantee sufficient 
degrees of freedom for permutation whereas at the same time 
not forfeiting the low frequency detection for the short data sets, 
which gets the more difficult the shorter the subset becomes. In 
order to facilitate proper summation of the results, the transfor- 
mations of the n subsets need to have the same frequency resolu- 
tion and frequency window so one really sums up the function at 
equal points. If the resulting periodogram still shows appreciably 
high;^'^ values, it can be deduced that the corresponding periods 
are actually physically relevant. If otherwise the peaks become 
increasingly broadened, it is rather likely that they do not corre- 
spond to a physically relevant period but stem from noise and/or 
aliasing effects. 

Another method is "Random-noise-modelling": We take the 
data points x,- and their error-bars 5,, create a uniformly dis- 
tributed random number e [-5,, -1-5,] and add ^, to the cor- 
responding Xi. The idea behind this is to weigh the influence of 
a bad signal-to-noise ratio. If a data set has good signal-to-noise 
ratio, the random-noise-modelling will scarcely effect the result- 
ing periodogram at all as the interval +5, is small and therefore 
the new set, created by the Random-noise-modelling is almost 
identical to the original one. If, however, the 6, are large, i.e. 
the data has a bad signal-to-noise ratio, the new set will dif- 
fer significantly from the old and it is to be expected that the 
period-analysis will yield different results. So we have found a 
way of assigning tendentially larger uncertainties to the results 
from data sets of bad quality . It might be argued that the dis- 
tribution of the ^, should be Gaussian, other than uniform, how- 
ever, we preferred the uniform distribution as it leads to larger 
differences of the new sets and the originals and gives therefore 
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a larger uncertainty in the final results, which we are otherwise 
underestimating. The variance of the new set is statistically con- 
verging towards that of the original as for each data-point x, the 
6i and from that the ^, is calculated. The distribution of the vari- 
ance is given in Fig.|6] This graph shows that the distribution of 
the variance, and thus of the rms of the simulated lightcurves, 
is gaussianly distributed around the original variance. We fitted 
with the formula 



1 



y '■ 



X exp 



10 



(InAl) 



(13) 



where fi - 0.306117 is the mean of the original set and - 
0.030765 a constant for normalisation. The fit yielded Aq - 
0.009955 with a correlation coefficient of 0.975010. Thus we 
are left with a new lightcurve, that will, statistically, have simi- 
lar rms and variance as the original one. We, again, pick n of the 
new lightcurves and execute the transformation. 

Throughout the analysis we took n = 3 representative sam- 
ples, each. 

The results are accumulated (= added up) onto one final pe- 
riodogram, which provides the means of determining the signif- 
icance level cr of the whole analysis. If now the periodogram 
resulting from the summation should exhibit its highest peak 
at places where the underlying single curves show nothing but 
a forest of little noise peaks, we deduce that the correspond- 
ing period is merely due to noise and has no physical meaning 
whatsoever This eff'ect is mainly observed for very small peri- 
ods {P < 10 days), whereas for higher periods the displacement 
of the large peaks along the period-axis dominate. 

Both random-methods were applied to all data sets indepen- 
dently, i.e. the original sets were taken, for each of them we cre- 
ated 3 sets using the Random-noise-modelling and (again from 
the original) 3 subsets using the random-sampling. The analy- 
sis of the data and its random models was carried out respect- 
ing all three numerical methods (see above). The simulated light 
curves are treated on equal footing as the original set, no weigh- 
ing is applied, nor are the curves mixed. Only the very final pe- 
riodograms (Fig. |7]i present the sum of all 7, individually anal- 
ysed, lightcurves for each experiment. As shown in Fig. [8] the 
Lomb-technique for high oversampling parameters may display 
oscillations, which appear to be Gaussian-shaped and more or 
less symmetrical around some central period. We believe these 
oscillations to be a mere artifact due to choosing the oversam- 
pling parameter slightly too big, a choice which is however nec- 
essary in order to resolve larger periods at all. Therefore we 
chose to fit enveloping Gaussian curves onto these structures. 
The cr levels were extracted from the HWHM value of the enve- 
lope, that was fitted onto the final (accumulated and normalised) 
periodograms. 

For the SigSpec tool, however, this formula is not applica- 
ble and also the method of fitting the envelope is generally un- 
derestimating the uncertainty. According to Kataoka| ( 2008) 1 the 
normalised power spectral density exhibits a power-law P(f) oc 
y-i.o fQj. (pink) frequencies, and the transformations were 
consequently corrected thus ("pink corrected"). In order to avoid 
spurious signatures from gaps, the data sets were each examined 
for their sampling structure so as to treat peaks arising around 
the positions of recurring data-gaps with additional caution. 

5. Results and Discussion 

In this section we present the results of our analysis methods 
applied to the various data sets available for Markarian 501. 



Distribution of Variance for Random-Noise-Modelling 

Histogram of 300 samples in 16 bins 
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Figure 6. Histogram of the variance of 300 Random-noise- 
modelling samples and the one-parameter fit onto a Gaussian- 
Distribution yielding a correlation coefficient of 0.975010 



First of all we show the periodograms calculated from the data 
(Figs .[7][8]l, th en we present the results of the hypothesis tests (in 



Figs. |i if- 14 1 and summarise all this in Table[T] 



In Fig. |7j and Fig. |8] we show the typical signatures of 
the fundamental frequency in RXTE and the first harmonic in 
WHIPPLE, respectively. The periodogram of the preliminary 
WHIPPLE data was produced using a very high oversampling 
parameter, so as to resolve the long periods; however the data 
set is too short for the Lomb-Scargle algorithm to produce re- 
liable results. In addition the Epoch-Folding-phase diagram in 
Fig. |9] shows the signature of the 72.6 day-period to be clearly 
sinusoidal. The overall hypothesis tests for both the first and sec- 
ond harmonic of 72 + 4.3 days are shown in Fig. 11-14 The 
hypothesis-tests are evaluated according to standard statistics us- 
ing two-sided tests on the specified confidence level. The points 
with their error-bars are the identified periods with their 1-cr- 
uncertainties whereas the boxes show the area in which the hy- 
pothesis can be accepted, if overlapping. We find, that the 23-day 
period can be accepted on a 5%-level for the SigSpec results 
(cf. ^ 



(cf. 



Fig. 
Fig. 



12 



11 



and on a 1%-level for the Lomb-Scargle results 
when all four experiments, that we consider sig- 
nificant, are compared. That means, that we accept a hypothe- 
sis even if the SWIFT results indicate the contrary (see Sec. |3] 
for discussion). For the Lomb-Scargle analysis of the 36-day 
period only the data-sets of sufficient length could be consid- 
ered due to the fact that the detection of low frequencies re- 
quires longer time spans in the original data. Fig. [13] shows that 
a period of 1/2 ■ (72.2 ± 2.3) days can be considered present 
in the two data sets on a 5%-level. In Fig. [14] we show that 
on a 5%-level, the first harmonic, 36 days, as identified by the 
SigSpec algorithm is still consistent within the 1-cr uncertainties 
of the fundamental period, 1/2 ■ (72.6 + 4.3) days, even though 
we probably underestimated the uncertainties in our method for 
this algorithm. The phase diagram in Fig. [9] shows that of the 



4-parameter fit onto Eq. 12 yields 0.621577 for the correlation 



coefficient, with an RMS per cent error of 0.18285 for the pa- 
rameter A2 = 0.0865454, which corresponds to a period of 72.6 
days via the fundamental relation oj - 2n/T. We interpret this 
high correlation coefficient as a clear sign the periodicity itself 
to be of sinusoidal shape. 

The approach to look for the harmonics of 72 days is moti- 
vated chiefly by Figs. [T] and [9] for this period dominates over all 
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RXTE 50901-54803: LombScargle Periodogram 

pink corrected 




20 30 40 50 60 70 80 90 100 HO 120 130 
Period [days] 



Phasediagram for Period = 72.6 

Epoch-Folding without any coiTcctions 




Phase [degree] 



Figure?. Pink-corrected Lomb-Scargle Periodogram for the 
RXTE data: 1998-2008 of Mrk 501. This final periodogi-am is 
the sum of the periodograms of the original lightcurve and the 
simulated lightcurves. The dashed graph is the periodogram of 
the original data set 

Whipple 06-08: Lomb-Scaigle Periodogram 

pink corrected 




Figures. Pink-corrected Lomb-Scargle Periodogram for the 
preliminai-y WHIPPLE data 2006-2008 of Mrk 501. This final 
periodogram was produced using a very high oversampling pa- 
rameter so as to resolve long periods in this short data set. Thus 
the influence of artifacts dominates for periods larger than 30 
days 



otiier structures. However, the ASM data is the only set, plentiful 
and long enough so the algorithms can detect and identify a pe- 
riod clearly. Given the fact that all other samples are dominated 
either by noise or gaps it is of great importance to seek for the 
persistent appearance of the harmonics throughout all sets, even 
though a single detection of a period in one sole set would for 
itself scarcely be considered significant. 



This study allows for a new limit on black hole masses in the 
possible BBH system in Mrk 501. These limits are in accordance 
with measurements with the total black hole mass derived with 
other techniques. For Mrk 501 the observed data corresponds to 
/ = 8 (where / = S muxiv) / S miniv) is the observed maximum to 



Figure 9. Phase Diagram of the period - 72.6 day in the RXTE 
1998-2008 data: The Signature is clearly that of a sinusoidal sig- 
nal. 

The fit with the formula: y = Aq+Ai ■sin(A2-xH-A3) yielded Aq = 
0.306917, Al = 0.051881, Aj = 0.086545, A3 = -0.466277 




log (m/lO^Mo) 

Figure 10. Required mass dependence for a BBH system in Mrk 
501. The dotted and dot-dashed lines show the upper and lower 
limits of the binary mass estimation derived from the M, - cr 
relation ( Woo et al.,2005 ), respectively. The solid (y^ - 15) and 
dashed (y^ = 10) thick curves are given by the Doppler condition 
for inclination angles / - lljb and an observed period of 72 
days. The dashed thin curves represent the same relation for a 
period of 23 days. The allowed mass-range has been indicated 
by filled areas. A TeV spectral index of a = 1.2 has been apphed 
for the calculation. 



minimum ratio of the spectral flux density) and a TeV spectral 
index of or ^ (1-2 - 1.7) (cf. Aharonian et al. [1999). Together 
with our proposed periodicity of 72 days we can determine the 
masses of the constituting black holes of the binary system by 
simultaneously solving Eq. ([8]) and (j9]l. The results are shown in 
Table [2] In Fig.[TO]we show the limiting black hole masses for 
the large as well for the small black hole. These are compared 
for the new found 72-day period and the previously found 23- 
day period. 

The obtained value for the separation of the BBH obviously 
changes in time as a consequence of gravitational radiation. In 
Fig. 10 we show the limiting black hole masses for the large as 



well for the small black hole. These are compared for the new 
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SWIFT 
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level 


Period [days] 


level 


Period [days] 


level 
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0.14 














36.9 ±2.3 


0.073 


32.7 ± 0.5 


0.029 


36.6 ±4.1 
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21.8 ± 1.2 


0.039 


21.6 ±0.3 


0.019 


22.4 ± 3.5 


1.5 


21.1 ±3.5 
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Table 1. The 5%- significance levels of the corresponding hypothesis tests brought face to face with the identified 72 day period 
and its higher harmonics 



Comparison on the 5% -Significance Level 

tile Hypothesis of 23 day periodicity can be accepted 



29 - 

28 - 



LombScargle with 1 -Sigma 
5% Significance Level 



Comparison on the 5% -Significance Level 

the Hypothesis of 23 day periodicity can be accepted if SWIFT is neglected 



26 - 
25 - 

24 - 

O 

S 21- 
D- 

20 - 
19 - 

18 - 

17 — 



SigSpec with 1 -Sigma 
5% Significance Level 



Figure 11. The two sided-hypothesis test of whether the 23- 
day period can be accepted within the 5%-level for the Lomb- 
Scargle algorithm: All found periodicities are consistent on the 
5% confidence level. 



Figure 12. The two sided-hypothesis test of the 23-day period 
for the SigSpec algorithm: The SWIFT data does not support the 
hypothesis, however, we very much doubt the quality of that data 
set and consider the hypothesis as valid on the 5%-level. 



found 72-day period and the previously found 23-day period. 



= 1/76 


1/10 


1/15 


mho'* Mo 




1.68 


0.52 


m[io'*Mo 




2.49 


3.65 


£/[l0'«^cm 




7.98 


13.70 


Pk [yrs] 




19.11 


42.96 



Table 2. The maximum masses of the constituting black holes, 
the separation d and the intrinsic orbital period for the inclination 
angles / ^ lljh w ith bulk Lorentz factor ji,. A total mass of 
Kf -^^MQ ( |Woo et aL.2005^ and a TeV spectral index of a = 1.2 
have been applied for the calculations. 



Our analysis shows two main results: Even for the quies- 
cent state of the source MrkSOl we find that on the one hand 
the well-established 23-day period (| Osone|[200 6 ) is still unde- 
feated, while on the other a 72-day period is upcoming. These 
results are not in contradiction since the 23- and also the 36-day 
periods can, within the uncertainties, be regarded as the second 
and first harmonic of the 72 ±4. 3-day period, respectively. Fig.|7] 
clearly identifies the 72 days, together with the persistent appear- 



ance of the first and second harmonics in the data sets (save for 
SWIFT, see discussion above). It is not fully clear how robust the 
applied methods are by themselves, however the combination of 
different approaches showing concordant results provides a solid 
argument in favour of the SMBBH model. 

This finding has possible consequences for our model of Mrk 
501: It allows for a new limit on the black hole masses in the 
source's possible BBH system. A limit which is firstly in ac- 
cordance with measurements of the total black hole mass as de- 
rived from other techniques and secondly allows for a larger ra- 
tio r - Mtotai/niBHi between the total mass of the binary and 
the mass of the companion BH. A scenario which is much more 



likely to exhibit the helical structure as proposed by Rieger & 
[Mannheim (2000 ). 

There is a reason why this has not yet been found: Short data 
sets are not able to show the 72-day feature since their length 
has to be well above these 72 days or aliasing will smear out this 
feature. 

Regarding the physical model there is now, with a multi- 
wavelength analysis, more reason to support the BBH model, 
but still long time observations, especially in very high energies, 
are needed to prove our results. 

The results of this study have an impact on x-ray, gamma- 
ray and gravitational wave astronomy: With evidence for multi- 
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Compai'ison on the 5% -Significance Level 

the Hypothesis of 36 day periodicity can be accepted 



-a 34 
o 



LombScargle with 1 -Sigma 
5% Significance Level 



Figure 13. The two sided-hypothesis test of the 36-day period 
for the Lomb-Scargle algorithm: the frequencies of the first har- 
monic are consistent on the 5%-level. The data sets from the 
experiments, which are not included, were too short for the al- 
gorithm to produce reliable results. 



Comparison on the 5% -Significance Level 

the Hypothesis of 36 day periodicity must be dicussed further 



O 
Oh 



X SigSpec with I -Sigma 
I 5% Significance Level 



Figure 14. The two sided-hypothesis test of 36-day period for 
the SigSpec algorithm: This period can be found to be persis- 
tent in the data on a 5%-level, if again SWIFT is not taken into 
account 



wavelength periodic variability, gravitational wave telescopes 
may help to tackle the problem of the origin of such variations 
and also aid in restricting the parameters of these sources. 
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